%%% File destined to calculate robust clustered standard errors for the
%%% main estimation files 


%%%%%%%%%%%% Risk parameters
fisher=hess_risk1j;


options = optimoptions('fminunc','FiniteDifferenceType','central', 'Algorithm','Quasi-Newton','Display','iter-detailed','TolFun',1e-5, 'TolX', 1e-7, 'HessUpdate', 'dfp', 'FiniteDifferenceStepSize',.1, 'MaxFunctionEvaluations', 100000, 'MaxIterations', 200);


coeff_vector=step1j_coeffs(size_factors+1:end);
n_para=size(coeff_vector,1);
%%% clustering is done at individual level so the calculations below are
%%% performed individual by individual
ind_total=1224;
ind_n=1;

temp_ind_id=ind_id;
temp_lottery_choice_ind=lottery_choice_ind;
temp_time_choice_ind=time_choice_ind;
temp_scores=scores;
temp_X_ind_s=X_ind_s;
temp_sim_factor_contribs=sim_factor_contribs;


%%% Records first derivatives for each individual. 
h=NaN(ind_total,n_para);
%%% records the precision of the derivative estimate
err_table=NaN(ind_total,n_para);

for ind=1:ind_total
    %%% Individual loop - get first derivatives for each individual with
    %%% respect to model parameters     
    ind_id=temp_ind_id((ind-1)*draw_n+1:(ind)*draw_n,:);
    lottery_choice_ind=temp_lottery_choice_ind((ind-1)*draw_n+1:(ind)*draw_n,:);
    time_choice_ind=temp_time_choice_ind((ind-1)*draw_n+1:(ind)*draw_n,:);
    scores=temp_scores((ind-1)*draw_n+1:(ind)*draw_n,:);
    X_ind_s=temp_X_ind_s((ind-1)*draw_n+1:(ind)*draw_n,:);
    sim_factor_contribs=temp_sim_factor_contribs((ind-1)*draw_n+1:(ind)*draw_n,:);
    
    [grad,err] = gradest(@ml_risk_score_types_cdf,coeff_vector);
    h(ind,:)=grad;
    err_table(ind,:)=err;
    ind
end

h_CRRA_cluster=h;

%%% N-Q adjustment 
B_CRRA_cluster=h_CRRA_cluster'*h_CRRA_cluster/(ind_total-n_para);
A=-fisher;
sandwich=inv(A)*B_CRRA_cluster*inv(A);
var_vector_robust_cluster=diag(sandwich)/ind_total;
se_vector_n_risk_robust_cluster=sqrt(var_vector_robust_cluster);
result_vector_risk_CRRA_robust_cluster=[coeff_vector se_vector_n_risk_robust_cluster];





%%%%%%%%%%%% Time parameters

coeff_vector_t=step2_coeffs;
n_para_t=size(coeff_vector_t,1);
fisher_t=hess_time;

ind_n=1;

temp_sim_theta_i=sim_theta_i;
temp_sim_kappa_i=sim_kappa_i;
temp_sim_1j_contribs=sim_1j_contribs;

%%% Records first derivatives for each individual. 
h_t=NaN(ind_total,n_para_t);
%%% records the precision of the derivative estimate
err_table_t=NaN(ind_total,n_para_t);

for ind=1:ind_total
    %%% Individual loop - get first derivatives for each individual with
    %%% respect to model parameters    
    ind_id=temp_ind_id((ind-1)*draw_n+1:(ind)*draw_n,:);
    lottery_choice_ind=temp_lottery_choice_ind((ind-1)*draw_n+1:(ind)*draw_n,:);
    time_choice_ind=temp_time_choice_ind((ind-1)*draw_n+1:(ind)*draw_n,:);
    scores=temp_scores((ind-1)*draw_n+1:(ind)*draw_n,:);
    X_ind_s=temp_X_ind_s((ind-1)*draw_n+1:(ind)*draw_n,:);
    sim_factor_contribs=temp_sim_factor_contribs((ind-1)*draw_n+1:(ind)*draw_n,:);
    
    sim_theta_i=temp_sim_theta_i((ind-1)*draw_n+1:(ind)*draw_n,:);
    sim_kappa_i=temp_sim_kappa_i((ind-1)*draw_n+1:(ind)*draw_n,:);
    sim_1j_contribs=temp_sim_1j_contribs((ind-1)*draw_n+1:(ind)*draw_n,:);
    
    [grad,err] = gradest(@ml_time_score_types_cdf,coeff_vector_t);
    h_t(ind,:)=grad;
    err_table_t(ind,:)=err;
    ind
end

h_CRRA_cluster=h_t;

%%% N-Q adjustment
B_CRRA_cluster_t=h_CRRA_cluster'*h_CRRA_cluster/(ind_total-n_para_t);
A=-fisher_t;
sandwich=inv(A)*B_CRRA_cluster_t*inv(A);
var_vector_robust_cluster_t=diag(sandwich)/ind_total;
se_vector_n_time_robust_cluster=sqrt(var_vector_robust_cluster_t);
result_vector_time_CRRA_robust_cluster=[coeff_vector_t se_vector_n_time_robust_cluster];


